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Abstract 

In this paper we present a new methodology for option pricing. 
The main idea consists to represent a generic probability distribution 
function (PDF) via a perturbative expansion around a given, simpler, 
PDF (typically a gaussian function) by matching moments of increas- 
ing order. Because, as shown in literature, the pricing of path depen- 
dent European options can be often reduced to recursive (or nested) 
one-dimensional integral calculations, the above perturbative moment 
expansion (PME) leads very quickly to excellent numerical solutions. 
In this paper, we present the basic ideas of the method and the relative 
applications to a variety of contracts, mainly: asian, reverse cliquet and 
barrier options. A comparison with other numerical techniques is also 
presented. 

Keywords: Option pricing, Barrier options, Asian options, Reverse Cli- 
quet options, Discrete Monitoring, Quadrature, Gram-Charlier Series, Black- 
Scholes. 

1 Introduction 

A central problem in today finance is to price correctly and rapidly exotic 
options. Indeed, in the last years, exotic options have became quite pop- 
ular and financial institutions traded larger and larger quantity of these 
sophisticated financial instruments. Moreover, it has became quite common 
to encounter exotic options embedded in structured bonds, which are ad- 
dressed also to small investors. There are no limits to market fantasy in 
the construction of new option varieties. It is therefore of great importance 
to develop secure and fast methodologies for option pricing. The literature 



on this argument is quite huge. Many approaches have been proposed, the 
most important of which are summarized below: 

(1) analytical solutions: usually they are limited to some special cases 
(simple contracts where the underlying asset is usually assumed to 
be governed by a log-normal process) j2J. Furthermore these exact 
solutions cannot be extended to more complex contracts. 

(2) binomial or trinomial tree: these popular approaches are exten- 
sively used both in academic literature [U |5j as well as financial in- 
stitutions because they can be easily adapted to many different situa- 
tions and are simple from a conceptual point of view. However many 
drawbacks are present: 1) the convergence is generally slow and the 
algorithm becomes rapidly inefficient in high dimension (e.g. for bar- 
rier options, the problem dimension scales linearly with the number 
of monitoring dates); 2) the error is not monotonic by increasing the 
number of nodes [Zj ; 3) usually the underlying movements are assumed 
to be described by a log-normal process; 4) the result accuracy depends 
on how the lattice is defined and last but not least, the method be- 
comes critical for particular choice of market data (e.g. for barrier 
options, when the spot price is close to barrier level 0). 

(3) Monte Carlo methods: MC is a powerful method for numerical 
pricing calculation jlOl I12j . Indeed when the problem dimension 
becomes high, it is the natural choice respect to binomial or trino- 
mial tree. That because the pricing error, in MC, scales as 1/V~N 
(where N is the number of simulations) independently from prob- 
lem dimension. Moreover, Monte Carlo can be implemented quite 
straightforward for all path dependent european options, it therefore 
not surprising that this method is largely used in banks and financial 
institutions. On the other hand, the convergence speed is low, mak- 
ing MC approach computationally demanding; that have originated 
a plethora of improved convergence methods and their relatives (re- 
duced variance techniques, quasi Monte Carlo approach based on low 
discrepancy sequences etc. [T2*] ^ 

(4) PDE solution: as well known, the problem of option pricing can be 
formulated in terms of partial differential equation (PDE) 2]; therefore 
a possible strategy, to price derivatives, consists to solve numerically a 
PDE. A good example of this approach, in the case of path dependent 
option with discrete sampling, is given in [1311141 ITS] . A strong point in 
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favor of this method, is the possibility of considering different processes 
alternative to geometric brownian motion (GBM). However some care 
must be taken in the implementation of the method. 

(5) quadrature method: the option pricing of path dependent european 
options could be reformulated as a path integral. In the case of discrete 
fixing dates, this path integral reduces to a multi dimensional integral 
(whose dimension correspond to the number of observation dates). As 
shown by some authors [El H3 HH H3 H3 [2T], often this multiple 
integral can be further reduced to a series of recursive (nested) one- 
dimensional integrals (one for each fixing dates), whose solution leads 
to an estimation of the option price. 

The problem faced by quadrature methods is how to treat accurately 
the density function involved in each of these one-dimensional nested 
integrals. Three main different approaches have been proposed in lit- 
erature: 

(i) a basic simple idea ^3 E] consists to discretize, via a grid of 
points, the density at each time step. Then, by adopting a re- 
cursive algorithm, the density function at the next time step is 
numerically calculated basing on the values at the previous step. 
This numerical scheme can prove to be expensive, from a com- 
putational point of view, as the result accuracy depends directly 
from the number of points composing the grid. 

(ii) alternative to crude Numerical Recursive Integration (NRI) pre- 
sented above, one can resort to parametric approximations of the 
real density, by finding the distribution parameters via quadratic 
fits. For instance, in the case of asian options, a popular choice 
is to parametrize the density of underlying arithmetic average, 
at each observation date, as a log-normal distribution. An im- 
provement was suggested by Lim |19| . who proposed a more gen- 
eral mixed density approximation in order to take in account the 
skewness of the distribution which comes up when high value of 
underlying volatility are considered. 

However the results, depend strictly on the arbitrary choice of a 
particular form used to model real densities. In other word, we 
do not have a rational criteria in order to individuate the better 
parametric representation. 

(iii) a third alternative consists to combine numerical integration with 
function approximation. In JH]) the authors consider the prob- 
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lem of a barrier option with discrete monitoring dates. They ap- 
proximate the solution at each observation date by a Chebyshev 
polynomials and solve the integral via Gauss-Legendre quadra- 
ture. An alternative approach has been proposed by Fusai et al. 
in again for barrier options. The authors approximate the 
solution thorough a linear combination of hat functions. All the 
task consists to relate the coefficients of the linear combination 
at time ti + \ to the coefficients previously computed at time U. 
They also provide an error estimate, which depends on the grid 
spacing used to interpolate the density function. This error is 
found to be quadratic in spatial discretization. 

In this paper we propose a new methodology which belongs to the area 
of quadrature algorithms. The method, we have called perturbative moment 
expansion (PME), focuses the attention on distribution moments instead of 
PDF's. The main idea is based on two legs: 

(a) a generic PDF can be always reconstructed, knowing its moments, by 
resorting to an extension of Gram-Charlier Series |22j . More precisely 
a PDF can be expressed as a perturbative series expansion around an- 
other (simpler) PDF, by matching all moments (up to a given order) of 
the original distribution. These arguments are presented in section El 

(b) on the other hand, basic operations involving PDF's, are simpler in 
terms of moments. For instance, in section |31 we show how the con- 
volution product between two distributions, can be reduced to simple 
arithmetics by reformulating the problem through moments. 

Because often the option pricing can be reduced to a recursive numerical 
integrations over density functions, PME permits to solve efficiently the 
problem with little computational efforts (section Indeed by including in 
the calculation just the first four moments (mean, volatility, skewness and 
kurtosis) the option pricing accuracy could be already high. 
In particular in section^we show how to implement our method for a variety 
of contingent claims, that is: barrier, asian and reverse cliquet options. 
Finally section |SJ is devoted to discuss briefly some conclusions and remarks. 
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2 Perturbative moments expansion of a probabil- 
ity distribution function 



Among quadrature methods for option pricing, a key problem is to model 
accurately PDF. In this section we will show how it is possible to accomplish 
the task in a very efficient way. 

Let us consider a generic stochastic variable x with PDF given by P. 
The P moment of order k, is defined as: 

{[x- < x >p} h ) P = / (x- < x > P ) k P(x)dx . (1) 

J — oo 

where < x > is the mean of P: {x)p = J x P(x) dx. 
The normalized moment of order k, is defined as follows: 

4" = llZ - - X >P]t)P . ■ CO 
{{[x- < X >p] 2 )p} 2 

By indicating with op the square root of second moment (i.e. the standard 
deviation of P): ap = {([x— < x >p] 2 )p}z, we can always decompose x as: 

x = < x >p +apy , 

m = J-p(^L^), (3) 
ap \ ap J 

where y is a stochastic variable with zero mean and unit variance and P 
denotes its PDF. 

It turns out that the probability density function P, can be represented as a 
series expansion in terms of a polynomial multiplied by the normal density 
$0,1 (with unit variance and zero mean). The polynomial accounts, in this 
way, the departure of the original PDF from normality. This expansion is 
known in literature as Gram-Charlier Series Usually the infinite series 
is truncated to a given order, very often up to fourth moment, while higher 
moment corrections are neglected. In such a way it is possible to incor- 
porate adjustments in the probability distribution for non-normal skewness 
and kurtosis effects (the former being the third moment and accounts for 
asymmetric tails, while the later corresponds to the fourth moment and in- 
corporates, for value higher than three, fatness in the tails). 
In past years some authors HHHiniUZl have resorted to Gram-Charlier Series 
to derive corrections to Black & Scholes formula for plain vanilla options, 
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by including skewness and leptokurticity effects in the distribution of stock 
returns. 

The Gram-Charlier Series, truncated to the fourth moment, reads: 



P(x) 



1 + f H 3 (x) + ^ H 4 (x) + 



*o,i(«) , (4) 



where $0,1 is the normal distribution with zero mean and unit variance, /U3 
and /X4 are respectively the skewness and kurtosis of P and {H n (x)} ^ 
denote the Hermite polynomials pQ. 
More generally, we can write: 



K 



P(x)« [YtCi J ) *o,i(x) , (5) 



u=0 



where the coefficients Cj (up to order if) can be easily computed by imposing 
the equivalence of the first K + 1 moments of both sides of equation ([5]>: 



$><^ +J W = (^>P > for: i = 0,l,2, K, (6) 

where we have indicated with < x l >$ 1 the gaussian moment of order i. 
Remember that for a gaussian distribution with zero mean and unit variance: 

< X s >$ x = if j is odd , 

< x 3 >$ j = 1 • 3 • 5 • • (j — 1) if j is even , (7) 

moreover, P is a PDF with zero mean and unit variance, therefore: 
< x° >p= 1, < x 1 >p= and < x 2 >p= 1. 

By reverting the linear system equations ©, the coefficients {cj}q' can be 
easily found in terms of a linear combination of the first K + 1 moments 
of P. Therefore, in some sense, the equation can be regarded as a 
perturbative moment expansion (PME) of P, around gaussian function $0,1 • 
The convergence of the series to P is guaranteed by some theorems |23| . 
It is straightforward to extend the equations (J5J) and ©, to the case where 
$o,i) is substituted with a generic PDF, ip (again with zero mean and unit 
variance). Of course, it is intuitively, that the quality of approximation 
in ©, at a given order K, depends on the "distance" between P and ip. 



6 



3 The convolution of probability distributions via 
moments calculation 



Let us consider two independent stochastic variables x\ and x 2 , each of 
theme characterized by a different probability distribution function P\ and 
P 2 . The sum of the two variables: z = x\ + x 2 is again a stochastic variable 
with a probability distribution, P, given by the convolution of P\ and P 2 : 



r r r+oo 

P(z) = / / Pi(x\) P 2 {x 2 ) dxi dx 2 = / Pi{xi) P 2 {z - xi) dxi . 

J JXl+X2=Z J —oo 

(8) 

Starting from the above equation, it easy to find the composition law con- 
necting the P moments to Pi and P 2 moments: 

(z) P = (x 1 ) Pl + (x 2 )p 2 




On the opposite of eq. the moments composition law is quite simple and 
straightforward. Therefore an alternative approach to compute the convolu- 
tion product between Pi and P 2 , in order to find out P, consists to compute 
P moments (up to a certain order) by using Q and then reverting to P by 
considering the perturbative expansion (J5J). The potentiality of this tech- 
nique could be appreciated by considering the following problem: given a 
stochastic variable z characterized by a probability distribution P, we ask 
to find out the PDF, P, such that the convolution: P\/ 2 o P\j 2 is equal to P 
(in other word we are ask to extract the square root of P respect to convo- 
lution product). This problem could be hard from a numerical point of view 
by considering the eq. (jHJ but becomes quite simple by resorting to PME 
and moments composition law. Indeed, within PME scheme, the problem 
reduces to simple algebra: once we have got the moments characterizing 
the square root of P, via inversion of eq. ©, it is simple to reconstruct the 
corresponding PDF by using equation ©. 

If we consider the sum of n i.i.d. variables, with distribution function 
Pi, the PDF of the sum, is given by P n = Yl\ Pi (where JT refers to the 
convolution product). In terms of moments this equation becomes: 

([x- < x >p n ] k ) Pn =n([x- < x > Pl ] k ) Pl + 
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x— < X >p 



k—h 



) P . ) ([x- < x > Pl ] u ) Pl . (10) 



Equation permits to evaluate iteratively, all P n moments. More pre- 
cisely, a solution can be easily found for k = 2, then, starting from the 
knowledge of k — 1-moment of P n , we have a rule to construct the next 
£;-moment of P n . 

As an example, we can solve iteratively equation (|10j) . for first values of k 
(i.e. k = 2, k = 3 and k = 4, which give the evolution laws, respectively, for 
volatility, skewness and kurtosis): 



{[x- < x > Pn ] ) Pn = n {[x- < x > Pl ) ) Pl , 



(3) 



u (4) 



(3) 
1/2 



r? 



3 + 



(4) 



(11) 
(12) 

(13) 



More generally, starting from eq. (|1U|1 . it possible to derive a perturbative 
series expansion of reduced moments [ip^ , in powers of . In equation (|14|) , 
we show the first corrections beyond gaussian result (n = oo): 



(k) 



(fc) 



where /i^ represents the normalized gaussian moment of order k. 
Formula (|14[l. shows that for large n, the PDF, P„, converges (according to 
central limit theorem) to a gaussian distribution and the first two corrections 
to the asymptotic limit, depends only from fjp^ and y£p — 3 (this result is the 
PME version of a famous theorem obtained by Kolmogorov and Gnedenko, 
for probability distributions, in 1954 [H] 1 ). 

lr The theorem provide an asymptotic series expansions of distribution function P n as: 



(fc) 



k 



, (3) 

1 Hp[ 



3! nV2 



+ o 



,( 4 ) 



3/2 



n 



I , Hk- 2)^-3 Q 



4! 



71 



k odd , 



k even , 



(14) 



Pn(x) - $o,i (*) = *o,i(a;) ^ 



7(3) 



(15) 



J=0 



where Qj is a polynomial whose coefficients depend uniquely from the first j + 2 moments 
of P n . The explicit form of {Qj} can be found in 12 II 
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4 Perturbative moment expansion for pricing Eu- 
ropean style options 



The option pricing of path dependent european options with discrete time 
monitoring, is equivalent to solve a multiple integral. As shown in litera- 
ture [HO El EDI Gil EH HU, often this multiple integral can be reduced to 
nested one-dimensional integrals (one for each fixing dates). The problem 
faced by quadrature methods is how to treat accurately the density func- 
tions involved in such integrals (see the introduction). 

In this paper, we propose a new scheme for modeling density functions, 
avoiding arbitrary choices of parametric functions as well as large grids which 
can lead easily to high computational costs. 

The basic idea relies from one hand on perturbative moment expansion of a 
generic PDF (equation ©) and from the other hand on moments composi- 
tion laws derived in section El 

We expose PME technique in practice, by considering three different kind 
of path dependent options: 

(i) asian options with discrete fixings (chapter 14. 

(ii) reverse cliquet options (chapter I4.2JI : 

(hi) barrier options with discrete monitoring dates (chapter 14.1 

for each of them, we present results obtained with PME approach. A com- 
parison with other popular techniques (Monte Carlo, Quasi Monte Carlo, 
recursive numerical integration, binomial tree etc.) is also reported. 



4.1 Asian options with discrete fixings 
4.1.1 Asian options: contract description 

We consider a call asian option written on a stock with initial value So an d 
volatility a. The interest rate, r, is assumed, for simplicity, to be a constant 
(indeed in our scheme this hypothesis is not necessary). 
The contract pay-off is defined as: 

Pay-off Asian = Max CHk^. - E , o) , (16) 
V m + 1 / 

where Si is the stock price at time t{ = i T/m, m is the number of equally 
spaced intervals and T indicates the time to maturity. 
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4.1.2 Asian options: PME problem formulation 

Equity price at time tj, can be written as: 

S j = S e^i=i( pi+Viei \ (17) 

where: e^s are independent iV(0, 1) random variables with zero mean and 
unit variance; pi and vf are respectively the mean and variance of stock 
log-returns calculated between and tf 

Pi = (*« — *»-i) > 

(18) 

Vi = a y/ti - ti-l . 

Now, we introduce the stochastic variables, ctj, defined recursively as follows: 

a = log(l + e ai ), (19) 
at = Pi + Vi ei + log (1 + e ai+1 ) , (20) 

Q>m — Pm ~\~ V m &m • (^^) 

aj's can be also re- written as 

a-i = Pi + ^ ii , (22) 

where pi and vf are the mean and variance of a« and {ii} are new stochastic 
variables with zero mean and unit variance. Ij's variables, on the opposite 
of ej's, are strictly related one to each other; more important, their PDF's, 
(referred in the following as PeJ, are in principle not normal. 
The option value, c asian , can be then calculated as: 



r+oo 

e~ rT / Max 

J— c 



e P0+v x 

m + 1 



E,0 



Pe (x)dx. (23) 



Therefore, the pricing of an asian option, with discrete fixings, is reduced 
to the calculation of a one-dimensional integral, with a PDF given by Pg . 
The problem is therefore reduced to compute recursively: pi, and P^ 
backward, until i = is reached. This task could be well accomplished, 
by using a perturbative moment expansion (PME) of P^ around normal 
distribution $0,1 ■ in detail, the iterative PME scheme for an asian option is 
defined as follows: 
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(A) fix the number of moments, I, to be used in PME algorithm; 

(B) let us start from i = m, Eqs. (|21|) and (|22|) give p m = p m , v m = v m 

and P £m = $0,1; 

(C) proceeding backward, for i < m, let us introduce a dummy stochastic 
variable: yi+\ = log ^1 + e^+ 1+ Vi+1 ei +^ . Knowing ?3j + i and 
P&i+n we can compute, numerically, all moments (up to order I) of 
Vi+i- 

Observing that aj is the sum of two independent stochastic variables 
(i.e. pi + Vi &i and j/j+i), we can calculate, by applying the moments 
composition rule derived in eq. ©: pi, Vi and all moments (up to I) 

Of Pe r 

Finally, by resorting to eq. ©, we are able to reconstruct P e . from its 
moments. 

(D) when i = 0, is reached, we can compute the option value via eq. (|23|). 

The higher the number of moments we include in the calculation, the higher 
will be the precision of option value estimate. However, with just including 
kurtosis and skewness (and neglecting higher moments corrections), we have 
got excellent results (the percentage error is less than 0.1 %, see tables ^ 
andHJ). 

4.1.3 Asian options: PME option pricing, a comparison with 
other techniques 

In order to compare our results with literature, we have considered two set 
of parameter values, reported in the captions of tabled and |2j 
In tabled we compare PME results (number of moments, I, ranging from 
4 up to 20) against other popular techniques: (i) Monte Carlo (MC) with 
iV = 10 8 scenarios, where we have resort to antithetic variables technique in 
order to reduce variance errors; (ii) Recursive Numerical Integration (RNI) 
and (iii) Mixed Density Approximation (MDA), which represents a modifica- 
tion of traditional RNI algorithm in which an appropriately parametrization 
of distribution functions is used. 

In table[21 a comparison of PME with MC, RNI and binomial tree (adopt- 
ing forward shooting grid algorithm) of Barraquand and Pudet [51 is also 
presented. 
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Asian option pricing with different techniques. 




17 = 5% 


<T = 10 % 


o — 30 % 


a = 50 % 


PME (Z = 4) 


4.30799 


4.90899 


8.79859 


12.96664 


PME (Z = 6) 


4.30798 


4.90899 


8.80149 


12.97995 


PME (Z = 8) 


4.30798 


4.90899 


8.80142 


12.98102 


PME (Z = 10) 


4.30798 


4.90899 


8.80153 


12.98124 


PME (Z = 20) 


4.30798 


4.90899 


8.80151 


12.98097 


MC 


4.30795 +/- 0.00003 


4.9089 +/- 0.00014 


8.80095 +/- 0.00062 


12.980 +/- 0.0012 


RNI 


4.308 


4.909 


8.801 


12.980 


MDA 


4.309 


4.911 


8.811 


12.979 



Table 1: Option pricing for an asian option with: So = 100, r = 9 %, T = 1 
year, E = 100 and m = 52 (number of intervals). The table shows option values 
computed with different techniques: PME - Perturbative Moment Expansion (with 
different values of I), MC - Monte Carlo (10 8 scenarios with antithetic variable 
technique), RNI - Recursive Numerical Integration (results are due to Lim |19p 
and MDA - Mixed Density Approximation (results are due to Lim |19| ) 

Finally, in figure ^ we show the relative error (in percent) between PME 
and Quasi Monte Carlo (QMC) results 2 . QMC estimates are based on 
N = 2 27 — 1 simulation series, making use of Sobol low-discrepancy se- 
quences. 

The results obtained, show that PME method give excellent results (error 
being less than 1%) also in the case we have retained few moments (I = 4, 
i.e. just kurtosis and skewness). 

The main advantage of PME respect to other numerical techniques are: 

(i) PME method, as all quadrature techniques, does not suffer from in- 
creasing fixing observations number or volatility values where, on the 

2 In spite of the fact that QMC method does not permit to estimate the error, we have 
preferred this algorithm to simple MC because a lower number of scenarios are required 
to obtain an estimate of a given precision. 
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Asian options pricing with different techniques. 




cr = 10 % 


a = 20 % 


<y = 40 % 


PME (( = 4) 


1.845289 


2.920988 


5.14583 


PME (! = 6) 


1.845298 


2.921097 


5.14677 


PME (( = 8) 


1.845297 


2.921095 


5.14677 


PME (i = 10) 


1.845298 


2.921097 


5.14678 


PME (i = 20) 


1.845298 


2.921097 


5.14678 


MC 


1.84517 +/- 0.00008 


2.9209 +/- 0.00018 


5.1463 +/- 0.00040 


RNI 


1.845 


2.921 


5.146 


FSG 


1.869 


2.960 


5.218 



Table 2: Option pricing for an asian option with: So = 100, r = 10 %, T = 
■J^r year, E = 100 and m = 91. The table shows option value computed with 
different approaches: PME - Perturbative Moment Expansion (with different values 
of I), MC - Monte Carlo (10 s simulations with antithetic variable technique), RNI 
- Recursive Numerical Integration (results are due to Lim ^H]), FSG - Forward 
Shooting Grid (Barraquand and Pudet in 0). 

opposite, the precision of MC estimates deteriorates when a or m in- 
crease. 

(ii) Once the final PDF is obtained, option prices can be easily computed 
for any So and strike price E, moreover that gives the opportunity to 
compute immediately the delta greek of the option. 

(iii) Respect to the methodology used in ^2], PME does not require to 
guess or make any hypothesis about the PDF form, nor to compute 
PDF on a grid of points. Indeed we have only to compute the first 
moments (e.g. for I = 4: mean, variance, kurtosis and skew, i.e. just 
four numbers) at each fixing date, making the option evaluation quite 
fast. For a general discussion about the advantage of PME algorithm, 
we remand to section [SJ 
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Err. PMEI4) 
Err. PME(6) 
Err. PME(8) 
Err. PME 10) 



fl fl & a 



-0.3 1 ' ' ' ' ' ' ' ' ' 1 

5 10 15 20 25 30 35 40 45 

Number of intervals 

Figure 1: Percentage error between QMC result (N = 2 27 - 1) and PME results 
(/ = 4-i- 10) for an asian option with: S = 100, r = 9 %, T = 1 year, E = 100, 
a = 50 % and m (fixing frequency) ranging from 4 up to 40. 



4.2 Reverse cliquet options 

4.2.1 Reverse cliquet options: contract description 

A reverse cliquet option is a product that typically have a maximum and 
minimum pay-out, which depends on the sum of negative performances of 
the underlying asset. 
A typical pay-out might be written as: 



Pay-off Rcv . d = Max 



(24) 



where the cliquet has m periods, Si is the underlying price at the end of the 
i'th period the strike level for the first period being SO (i.e. the 

spot price). L is the global floor (usually set to zero) and corresponds to 
the minimum amount the investor will receive at the expiration date. H is 
the maximum coupon payable by the option. 

The risk profile embedded in this option is clear: the investor is long the 
skew (if the skewness of stock price returns increases, it is less probable to 
have negative performances and therefore the contract value increases) and 
short volatility (indeed if volatility grows the contract price declines) . As we 
will see, the PME approach permits to calculate at each step the effective 
volatility and skewness of the sum of the underlying negative performances 
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(also in presence of non log-normal process), making more clear the risk 
profile. 



4.2.2 Reverse cliquet options: PME problem formulation 



the 



The problem formulation in terms of perturbative moment expansions, be 
comes in this particular case quite simple. 

Let us introduce the stochastic variables, Ri = Min ( ^'g,' 5 ''" 1 ; 0) , i 
minimum value between stock performance and zero. If we indicate with 
Pi the PDF of variables Ri, the density function of the sum of all negative 
performances (J2iL\ Ri), is given by: 



(25) 



i=l 



where the symbol, JJ, refers to convolution product. 

The equation (|25|). which is indeed a multiple integral, can be simplified, 
by reformulating the problem in terms of PDF's moments (see section EJ. 
More precisely, the P moments can be computed by combining, iteratively, 
the Pi moments via equation Q. On the other hand, the P{ moments can 
be easily calculated as follows: 



Ri(x) = Min 



su 



5ti x 



1,0 



x ~ N(0,1) 



(Ri) Pt = J Ri(x)$ ,i(x)dx, (26) 
([Ri- < R t >pf) Pt = J [Ri(x)- < Ri > Pi ] k $o,i(3 



(x) dx . 



Observe that the above integrals, can be easily "solved" in terms of a sum 
of cumulative normal functions. 

In equation ()26|) . $0,1 ( x ) represents the normal density distribution with 
zero mean and unit variance and 5ti = ti — ti—\. 

Once the P moments have been calculated, by resorting to eq. © , it is 
possible to reconstruct the corresponding probability distribution function 
as a series expansion around normal density, $0,1 • Then, the option price 
can be easily calculated as: 



option price Rov C1 
Note that: 



e" rT Max 



L, J (H + x) P(x) dx 



(27) 
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(i) if time intervals are constant the variables Ri are i.i.d. and therefore 
their moments can be computed once, making the option price evalua- 
tion quite fast (the laws governing moments evolution for i.i.d. permits 
to compute P moments by simple arithmetics, see equations (|10j) . (fl"2*|) 
and (|H1) . 

(ii) the inclusion of kurtosis and skewness effects in the process governing 
stock price variations does not require any modification to the algo- 
rithm, indeed it would be sufficient to substitute in eq. (|26j) to $0,1 
the PDF of a non log-normal process. 

(iii) by increasing the number of intervals, according to equation (|14jl. the 
PDF describing the sum of negative performances, converges to a gaus- 
sian distribution (see eq. dJ). We can expect therefore, that the 
Gram-Char lier series truncation, used in eq. ()27|) to model P(x), will 
become asymptotically exact for large value of m. 

4.2.3 Reverse cliquet options: PME option pricing, a comparison 
with other techniques 

In table we present PME results compared with QMC for different num- 
bers of equally spaced intervals, m, (ranging from 4 up to 36). The time 
intervals are maintained fixed to 1/12 year. In order to make results com- 
parable, we have considered H = mh {h constant) and option prices (in %) 
rescaled by a factor 1/m. 

PME results reported in table |3J show an excellent agreement with QMC 
estimates, at least for m > 12, even retaining only few moments in the per- 
turbative expansion (the error is less than 0.1 % with just the first four 
moments) . 

In figure 12 we plot relative error between PME results (I, ranging from 
6 up to 12) and QMC for different values of m, keeping the time interval 
constant. As we have stated in previous paragraph, the PDF of sum of 
negative performances converges, in the limit m — > 00, to a gaussian distri- 
bution, making PME asymptotically exact. As a consequence the agreement 
between PME and QMC becomes better and better increasing the number 
of intervals m. 
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Reverse cliquet options pricing 
with different techniques. 




m — 4 


m = 12 


m — 24 


m — 36 


PME (i = 4) 


1.41681 % 


1.01701 % 


0.82935 % 


0.72502 % 


PME (i = 6) 


1.42984 % 


1.01755 % 


0.82890 % 


0.72465 % 


PME (i = 8) 


1.43915 % 


1.01896 % 


0.82914 % 


0.72470 % 


PME (i = 10) 


1.43730 % 


1.01769 % 


0.82891 % 


0.72467 % 


PME (i = 20) 


1.43417 % 


1.01798 % 


0.82898 % 


0.72469 % 


QMC 


1.4336 % 


1.0179 % 


0.82898 % 


0.72469 % 



Table 3: Option pricing for a reverse cliquet option with: So = 100, r = 9 %, 
a = 30 %, 5 t = U- U-x = 1/12 year, L = and m = 4, 12, 24, 36. In order to 
make results comparable, we have considered H = rah (h — 4 %) and option prices 
(expressed in %) have been rescaled by a factor 1/m. The table shows option values 
computed with different techniques: PME - Perturbative Moment Expansion (with 
different values of I) and QMC - Quasi Monte Carlo (2 27 — 1 simulations). 

4.3 Barrier options with discrete monitoring dates 
4.3.1 Barrier options: contract description 

The pay-off of a generic knock out barrier option can be written as follows: 



Pay-off E 



_ J ^ 7 [5'(T)] if S(U) > Bi for all observation dates {U} 
1 b otherwise 

(28) 

where: J 7 is a generic positive function which represents the option pay-off if 
the barrier has not been touched, b is known as "rebate" and represents the 
amount paid by the contract if the barrier has been touched (often 6 = 0), 

Bi is the barrier level at time tj 3 , {tj}j=i,2, n are the discrete observation 

dates set, Sq the spot price, S(ti) the underlying price at time tj. 

For constant barrier level Bi — B Vi, in the present discussion we do not make any 
assumption about {Bi}. 
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S 0.2 



Err. PMEI6) 
Err. PME(8) 
Err. PME(10) 
Err. PME 12 



I - 5 I * • I i - i » » 



Number of intervals, m 



Figure 2: Reverse cliquet option: the graph shows the relative error between 
PME and QMC estimates, versus the number of intervals, m. Different values of 
I are considered (ranging from 6 up to 12). Data refers to: So = 100, r = 9 %, 
S t = ti — ti_i = 1/12 year, a = 30 %, L = 0, H = m h with h = 4 % and m = 
ranging from 4 up to to 40. 



4.3.2 Barrier options: PME problem formulation 

The problem implicit in a barrier option can be reformulated in terms of 
PDF's. The stock price a time ti is given by: 

Sfc) = S e^ 3=1 \ 2 ) 3 V 3 3 , (29) 

where 5tj = tj—tj-i and {ej} are independent normally distributed random 
variables with zero mean and unit variance. 
Let us define the following stochastic variables: 

Zl = (log ^ | S(tj) > Bj Vj < ij , (30) 

which embedding the condition that all observations until time ti are above 
the barrier level. 
Let us introduce: 

Pi = probability that: S(tj) > Bj Vj < i , (31) 
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i.e. the probability that the option is alive at time £j. 

Then, the option price can be expressed in terms of PDF of z m , P Zm , as: 



option price KO = e r( * m * o) 



[1-Pm)b + Pm J F(Soe x )Pz m (x)dx . 

(32) 

The stochastic variables {zi} can be recursively related one to each other, 
by the following equations: 



z = 0, 

Wi = z^i +17" -— I 5ti + a y/ 5ti ei , (33) 



a 2 



Zi = (wi | S(ti) > Bi) = (wi | Wi > log , 

where we have introduced the dummy variables Wi. 

The equations set (|33|) can be translated in terms of PDF's as follows: 

P Z0 (x) = 8(x), (34) 

Pwi = Pzi-! ° ®(r-a 2 /2)6ti,a 2 5t, , (35) 

P,(x) = { * P <«to if f> lo gf (36) 
lK ' \ otherwise v ' 

where: 5(x) is the Dirac delta distribution, "o" indicates the convolution 
product, < I ) p i „2 is the normal distribution with mean p and variance v 2 and 
hi is the probability to find Wi below log ^: 

("OO 

h= B P Wl {x)dx. (37) 

In this framework, the probabilities pi are then related by the recursive 
relation: 

Pi =p i - l -k i . (38) 

The PME iterative scheme for a barrier option evaluation can be therefore 
structured as follows: 

(A) fix the number of moments, I, to be used in PME algorithm; 

(B) repeat, starting from i = 0, the steps B1-B3 until i = m is reached: 
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(Bl) if i = 0, P zo is the Dirac delta distribution, therefore its moments, 
at any order are identically zero apart from moment of order zero 
which is, by definition, equal to 1. 

Otherwise for i > 0: knowing the moments of P Zi _ t , compute all 
moments, up to order I, of PDF P w . by using equation . 

(B2) by reverting to equations (JHJ), © and ©, reconstruct PDF P Wi 
by a Gram Charlier series around normal density function $o,i- 

(B3) given P Wi , from eq. (|3*B|) . compute all moments of P Zi up to order 
I. 

(C) Reached i = m, we know all P Zm moments up to I, again, by using 
equations ©, © and @ we can develop P Zm as a Gram Charlier 
series around normal density. One that is done, the option price can 
be easily evaluate by means of eq. (f32|) . 

It interestingly to note that within the scheme presented above, we need 
to evaluate only integrals of the form: 



which can be reduced to the evaluation of the well known inverse cumula- 
tive normal function 5 . Therefore we do not need any numerical methods to 
compute the required integrals. 

A last remark on PME scheme: if we want to treat also non log-normal pro- 
cesses for equity prices evolution, the only modification to be implemented 
regards the equations (|3*3*|) and (|35j) . where in place of ei and $o,i we must 
consider a non normal stochastic variable/PDF. Hence, when computing the 

4 The moments of a gaussian distribution are known at every order, see equation 0. 
5 For instance: 




(39) 



Io(d) 



Err_f(d) , 



h(d) 



1 _i 
i= e 2 



h{d) 



h{d) 



Err_f(d) - d 



(2 +d 2 ) 



1 
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P Wi moments via eq. (JSJ, we just need to consider the real PDF moments, 
different from those of a normal distribution. As an example in figure |3J we 
report how the option price changes when fat tails (i.e. non-normal kurto- 
sis) are considered. As expected, increasing leptokurticity of stock returns 



y 6.265 
J- 6.26 
6.255 
6.25 



Kurtosis of stock returns PDF 



Figure 3: Price of a down-out call option, calculated with PME (I = 20) by 
considering a leptokurtic behavior (i.e. kurtosis grater than 3) in probability density 
function of stock returns. So = 100, r = 10 %, T = 0.2 year, E = 100, a = 30 %, 
B = 89 and m = 5. 



distribution has the effect to enlarger the chances of touch the barrier before 
maturity, lowering therefore the option value. 

4.3.3 Barrier options: PME option pricing, a comparison with 
other techniques 

In order to make a quantitative comparison with other numerical methods, 
we have considered a down-out barrier option with constant barrier B. The 
pay-off at maturity is given by: 

p a Q fr _J S — E if S > E (40) 
1 otherwise 

i.e. the usual call plain vanilla pay-off. 

For this kind of option, recently Fusai et al. [H] have derived a closed form 
solution. Table 0] shows, for different values of barrier level, B, and obser- 
vation dates number, m, a comparison of PME with the exact result and 
other numerical techniques. 
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Barrier option pricing with different techniques. 


B 


89 


95 


97 


99 


89 


95 


97 


99 


m 


5 


5 


5 


5 


25 


25 


25 


25 


PME(4) 


6.27407 


5.68309 


5.17957 


4.49928 


6.20762 


5.11871 


4.14616 


2.82893 


PME(12) 


6.27954 


5.67106 


5.16749 


4.48947 


6.19971 


5.08100 


4.11596 


2.81157 


PME(20) 


6.28077 


5.67109 


5.16726 


4.48918 


6.20850 


5.08032 


4.11508 


2.81154 


PME(32) 


6.28076 


5.67110 


5.16725 


4.48917 


6.21032 


5.08096 


4.11559 


2.81224 


Exact 


6.28076 


5.67111 


5.16725 


4.48917 


6.20995 


5.08142 


4.11582 


2.81244 


QMC 


6.28075 


5.67111 


5.16726 


4.48912 


6.210005 


5.08156 


4.11561 


2.81233 


NRI 


6.2763 


5.6667 


5.1628 


4.4848 


6.2003 


5.0719 


4.1064 


2.8036 


CMF 


6.284 


5.646 


5.028 


4.050 


6.210 


5.084 


4.113 


2.673 


TT 


6.281 


5.671 


5.167 


4.489 


6.21 


5.081 


4.115 


2.812 


SRQM 


6.2809 


5.6712 


5.1675 


4.4894 


6.2101 


5.0815 


4.11598 


2.8128 



Table 4: Option pricing for a barrier option with: So = 100, E = 100, r = 10 
%, T = 0.2 year, a — 30 %. The table shows option price computed with different 
techniques: PME (I = 4, 12, 20 and 32), Exact results due to Fusai 0, QMC - 
Quasi Monte Carlo method (calculated with 2 27 — 1 scenarios), NRI - Numerical 
Recursive Integration (results are due to Aitsahlia 16 ), CMF - Continuous Moni- 
toring Formula (Broadie in |H]), TT - trinomial tree (results are due to Broadie 
and SRQM - Simpson Recursive Quadrature Method (Fusai in |21p. 

The agreement is quite good and it is better, in general, than the cor- 
responding result achievable by simple NRI, at least when the number of 
observation dates is relatively small (as it is in the example considered in 
table 

In figure 0J we report the relative percentage error between PME estimates 
and the exact result, for two different monitoring frequencies. The graph 
shows the accuracy of the PME method, which is already good (errors less 
than 0.1 %) just including few moments. 

However the number of moments to be considered in order to reach a good 
precision is higher respect to asian options (see chapter [4.1j) or reverse cliquet 
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5 10 15 20 25 30 35 

I - number of moments useO in PME 



Figure 4: Percentage error between PME estimates (I — 4, 6, 34) and exact 
result for a down-out call option with S = 100, r = 10 %, T = 0.2 year, E = 100, 
cr = 30 %, B = 99 and m = 5, 25. 

options (see chapter l4.2j) . The reason being the strong skewness (especially 
when equity spot price is near to the barrier or m is large) characteriz- 
ing the underlying probability distribution, P Zm , which leads to a strong 
departure from gaussian behavior. Indeed large skewness values lead to a 
Gram-Charlier expansions which may not be non-negative, especially on the 
distribution tails (compromising therefore the acceptance of the function as 
a true density). Because in a barrier option a key point is to measure the 
probability to touch the barrier at each step (i.e. to calculate accurately an 
integral on the left tail) the non positive definiteness of the PDF's tails, may 
lead to errors that compromise the pricing quality. That requires, especially 
when m is large, a correction to simple PME scheme, in order to fix (at 
least partially) the problem. A simple trick consists to substitute the step 
(B2) of the algorithm described in the previous paragraph, by a two layer 
procedure: 



(LI) reconstruct PDF, P w ., as a perturbative series expansion around 



(L2) (a) find the lower negative value xq such that: X)*=o c% Xq = 0; 



(B2) 



normal density: P, 




(b) considering a new function P Wi (x) 
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Due to the truncation procedure, we can expect that P Wi mo- 
ments depart slightly from initial target moments. However we 
can apply again the perturbative expansion procedure using as 
basic PDF the P Wi function 6 instead of a simple gaussian. 

In that way it is possible to derive an improved perturbative series ex- 
pansion, where the problem of non positiveness is partially fixed . 
In table [3 we report PME valuations (using such improved algorithm and 
including the first twelve moments, I = 12) versus exact results. The error is 
typically lower than 0.1 %, also for large number of observation dates, m. 



Barrier options pricing 

(improved PME / exact results) 


m 


PME(12) 


Exact result 


Relative error 


10 


4.1820 


4.18224 


0.006 % 


50 


3.1260 


3.12633 


0.011 % 


80 


2.9391 


2.93918 


0.003 % 


100 


2.8643 


2.86442 


0.004 % 


120 


2.8087 


2.80903 


0.012 % 


150 


2.7461 


2.7474 


0.047 % 


180 


2.7020 


2.70163 


-0.014 % 


200 


2.6746 


2.67682 


0.083 % 


220 


2.6531 


2.65545 


0.089 % 


250 


2.6253 


2.628099 


0.107 % 


280 


2.6018 


2.60534 


0.136 % 


300 


2.5879 


2.59056 


0.103 % 


500 


2.5018 


2.50259 


0.032 % 



Table 5: Option pricing for a barrier option with: spot price = 100, strike = 100, 
barrier = 98, r = 10 %, T = 0.2 year, a = 30 %. The table shows a comparison 
between the option value computed with an improved PME technique (I = 12) 
against the exact result (Fusai et al. in 0) for different m values. 



6 Before that, it is necessary to rescale Pwi, according to eq. JjJ, in order to obtain a 
PDF with zero mean and unit variance. 

7 Indeed, the above procedure does not guarantee that the new expansion is non neg- 
ative for any value of x. However numerical evidence shows a good improvement and a 
significantly reduction of the non positiveness problem. 
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5 Summary and Conclusions 



In this article we have proposed a new methodology for option pricing be- 
longing to quadrature techniques. Respect to other quadrature algorithms, 
where for instance density functions are modeled through a grid of points 
or polynomial interpolations, our method adapts a GramCharlier series ex- 
pansion around a given distribution (usually a gaussian function). 
The highlights features of the method are the following: 

(i) by representing PDF's through their moments, we are able to capture 
all the essential features of the problem by using just few parameters, 
improving computational performances and memory usage. How we 
have shown in the paper, the crude choice of retaining only the first 
four moments (i.e. kurtosis and skewness) in PME scheme, provides an 
excellent approximation of the option value (with errors less than 0.1% 
for asian and reverse cliquet options). Furthermore, by increasing the 
number of moments retained in PDF expansion, the numerical solution 
converges to the exact result, making possible, in principle, to increase 
progressively the estimate precision. 

(ii) The method is extremely simple and natural from a conceptual point 
of view and therefore easily implementable. 

(iii) It is extensible to different pay-off contracts (indeed the programming 
code needs very few modifications changing the contract features) . 

(iv) non log-normal processes (i.e. stochastic processes characterized by a 
non log-normal PDF), can be naturally treated within a PME scheme, 
without any modification to the algorithm. 

Moreover, the proposed method, maintains all the typical advantages of 
quadrature methods, that is: 

(i) we need to perform computations only at trigger times; 

(ii) the CPU time scales linearly with the number of observation dates. 

(iii) the price estimate is not too sensitive to changes in volatility; 

(iv) there are no time discretization errors. 

(v) it is easily to incorporate additional exoticity in the contract (e.g. for a 
barrier option, it is straightforward to include a time varying barrier). 
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